Compute interaction between Sleep deprivation effect and genotype effect at QTL peaks
Read data
File<-"../Data/Metabolite_BXD.txt"
x<-read.table(File,header=T,sep="\t")
n<-x[1,] # Contains the metabolites type
names(n)<-gsub(".","_",names(n),fixed = TRUE) # reformate a bit the data
x<-x[-1,] # We remove the first lines containing the metabolites type
Some of our results contains Null values (-1), we replace them by NA values
head(x[1:10,1:10])
x[x==-1]<-NA
# Take only BXD
x<-x[-which(x$Group %in% c("B6D2F1","C57BL6","D2B6F1","DBA")),]
Read FC-mQTL results
fcmQTL<-read.table("../../QTL/Metabolite/FC/QTLlist.txt",header=T)
Get marker
markers<-read.table("Genotype.FormatedName.geno",header=T)
n<-markers$Locus
colnames(markers)
## [1] "Chr" "Locus" "cM" "Mb" "BXD005" "BXD029" "BXD032"
## [8] "BXD043" "BXD044" "BXD045" "BXD048" "BXD098" "BXD049" "BXD050"
## [15] "BXD051" "BXD055" "BXD056" "BXD061" "BXD063" "BXD064" "BXD065"
## [22] "BXD097" "BXD066" "BXD067" "BXD070" "BXD071" "BXD073" "BXD103"
## [29] "BXD075" "BXD079" "BXD081" "BXD083" "BXD084" "BXD085" "BXD087"
## [36] "BXD089" "BXD090" "BXD095" "BXD096" "BXD100" "BXD101" "B61"
## [43] "DB1"
colnames(markers)<-gsub("BXD0","",colnames(markers))
colnames(markers)<-gsub("BXD","",colnames(markers))
colnames(markers)<-paste("BXD",colnames(markers),sep='')
markers<-markers[,colnames(markers) %in% unique(x$Group)]
# order marker
markers<-markers[,order(match(colnames(markers),as.character(unique(x$Group))))]
rownames(markers)<-n
Recompute quickly to get the peak
and then compute the interaction between genotype and sleep deprivation at QTL peak
library(qtl)
QTL<-read.cross("csv",file="../../QTL/Metabolite/FC/MetabolitesFC.QTLR",genotypes=c("BB", "BD", "DD","EE","FF"),allele=c("B", "D"),estimate.map=FALSE,na.strings=c("NA","-"))
## --Read the following data:
## 33 individuals
## 11186 markers
## 126 phenotypes
## Warning in fixXgeno.f2(cross, alleles): --Omitting 75 male heterozygote
## genotypes on the X chromosome.
## Warning in summary.cross(cross): Some markers at the same position on chr
## 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,X; use jittermap().
## --Cross type: f2
QTL<-convert2risib(QTL)
## Warning in convert2risib(QTL): Omitting 6411 genotypes with code==2.
QTL<-jittermap(QTL,amount=1e-6)
QTL <- calc.genoprob(QTL, step=1,error.prob=0.05)
intpvals<-c()
effects<-c()
topmarkers<-c()
QTLinfo<-matrix(nrow=0,ncol=10)
SavedMarker<-matrix(ncol=33,nrow=0)
colnames(SavedMarker)<-colnames(markers)
markers<-markers[,order(match(colnames(markers),as.character(unique(x$Group))))]
for (k in as.character(unique(fcmQTL$Phenotype))){
i<-which(fcmQTL$Phenotype==k & fcmQTL$Lodscore == max(fcmQTL[fcmQTL$Phenotype==k,"Lodscore"]))
id<-fcmQTL[i,"PhenotypeNumber"]
QTL.scan<-scanone(QTL,method="em",pheno.col=id)
st<-which(rownames(QTL.scan) == fcmQTL[i,"MarkerStart"])
ed<-which(rownames(QTL.scan) == fcmQTL[i,"MarkerEnd"])
QTLscan_small <- QTL.scan[st:ed,]
# Markers at QTL peak
markersQTL<-QTLscan_small[QTLscan_small$lod==max(QTLscan_small$lod),]
#remove pseudomarker
gr<-grep("c\\d+\\.loc\\d+",rownames(markersQTL),perl=T)
if (length(gr)>0){
markersQTL_Nopseudo<-rownames(markersQTL)[-gr]
}else{
markersQTL_Nopseudo<-rownames(markersQTL)
}
#QTLscan_small<-QTLscan_small[- grep("c\\d+\\.loc\\d+",markersQTL,perl=T),]
# if only pseudo marker available:
if (length(markersQTL_Nopseudo) == 0){
#estimate most probable genotype
marker<-rownames(markersQTL)[1]
Chrom<-fcmQTL[i,"QTL_Chromosome"]
genoproDD<-pull.genoprob(QTL, chr=Chrom)[,paste(marker,':DD',sep='')]
id<-as.character(QTL$pheno$Mouse_ID)
gm<-rep("B",length(id))
names(gm)<-id
gm[genoproDD>0.5]<-'D'
# A normal marker is available
}else {
marker <- markersQTL_Nopseudo[1]
gm<-as.character(as.matrix(markers[marker,]))
names(gm)<-colnames(markers)
}
if (nrow(QTLscan_small) > 0){
maxid<-which(QTLscan_small[,"lod"]==max(QTLscan_small[,"lod"]))[1]
marker<-rownames(QTLscan_small)[maxid]
gms<-gm[order(match(names(gm),colnames(SavedMarker)))]
SavedMarker<-rbind(SavedMarker,gms)
rownames(SavedMarker)[nrow(SavedMarker)]<-marker
geno<-gm[as.character(x$Group)]
geno<-as.factor(geno)
geno<-relevel(geno, ref="B")
Pheno<-fcmQTL[i,"Phenotype"]
Pheno<-gsub("_",'.',Pheno)
idxmeta<-which(colnames(x) == Pheno)
fit <-lm(as.numeric(as.character(x[,idxmeta]))~0+geno/factor(x$Group)+geno*factor(x$Condition))
Bestim<-summary(fit)[[4]]["genoB","Estimate"]
Destim<-summary(fit)[[4]]["genoD","Estimate"]
SDestim<-summary(fit)[[4]]["factor(x$Condition)SD","Estimate"]
DSDestim<-summary(fit)[[4]]["genoD:factor(x$Condition)SD","Estimate"]
intpval<-summary(fit)[[4]]["genoD:factor(x$Condition)SD","Pr(>|t|)"]
SD_B<-log2(Bestim+SDestim)-log2(Bestim)
SD_D<-log2(Destim+SDestim+DSDestim)-log2(Destim)
FCeff<-SD_B-SD_D
effects<-c(effects,FCeff)
intpvals<-c(intpvals,intpval)
topmarkers<-c(topmarkers,marker)
} else {
intpvals<-c(intpvals,1)
effects<-c(effects,0)
topmarkers<-c(topmarkers,NA)
}
QTLinfo<-rbind(QTLinfo,fcmQTL[i,])
}
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.
nrow(QTLinfo)
## [1] 89
QTL.scan[marker,]
intpvals<-p.adjust(intpvals,method="fdr")
QTLinfo$pvalint<-intpvals
QTLinfo$intlogFC<-effects
QTLinfo$topmarker<-topmarkers
fctest<-QTLinfo[,c(1,5,7,11,12,13)]
fctest[order(abs(fctest$intlogFC),decreasing = T),]
write.table(QTLinfo,file="Genotype_SD_Interaction_Metabolites.txt",sep="\t",quote=F,row.names=F)
plotFCMetabo<-function(WithVariant,WithoutVariant,Metabo,main){
File<-"../../Rdata/Metabolites_BXD_alldata.txt"
x<-read.table(File,header=T,sep="\t")
x[x==-1]<-NA
x<-x[-1,]
m<-x[,-c(1,2,5)]
m$Group
#WithVariant<-c("DBA",WithVariant)
#WithoutVariant<-c(WithoutVariant,"C57BL6")
MetaNSDMean<-matrix(ncol=1,nrow=0)
for (i in unique(m$Group)){
MetaNSDMean<-rbind(MetaNSDMean,mean(as.numeric(as.character(m[m$Group==i & m$Condition=="NSD",Metabo])),na.rm=T))
}
rownames(MetaNSDMean)<-unique(m$Group)
MetaSDMean<-matrix(ncol=1,nrow=0)
for (i in unique(m$Group)){
MetaSDMean<-rbind(MetaSDMean,mean(as.numeric(as.character(m[m$Group==i & m$Condition=="SD",Metabo])),na.rm=T))
}
rownames(MetaSDMean)<-unique(m$Group)
FC<-MetaSDMean/MetaNSDMean
d<-c()
col<-c()
geno<-c()
for (i in rownames(FC)){
if (i %in% WithoutVariant){
d<-c(d,FC[i,])
col<-c(col,rgb(0.5,0.5,0.5,0.8))
geno<-c(geno,"B")
}
if (i %in% WithVariant){
d<-c(d,FC[i,])
col<-c(col,rgb(0.85,0.67,0.44,0.8))
geno<-c(geno,"D")
}
if (i =="C57BL6"){
d<-c(d,FC[i,])
col<-c(col,rgb(0,0,0,1))
geno<-c(geno,"B")
}
if (i =="DBA"){
d<-c(d,FC[i,])
col<-c(col,rgb(0.7,0.5,0.3,1))
geno<-c(geno,"D")
}
if (i == 'D2B6F1' | i == 'B6D2F1'){
d<-c(d,FC[i,])
col<-c(col,rgb(1,1,1,1))
geno<-c(geno,"BD")
}
}
names(geno)<-names(d)
barplot(height=log2(d[order(d)]),las=1,horiz=TRUE,col=col[order(d)],main=main,xlab="log Fold-Change")
}
QTLinfo<-QTLinfo[order(abs(QTLinfo$intlogFC),decreasing = T),]
par(mfrow=c(2,2))
for (i in 1:nrow(QTLinfo)){
if (QTLinfo[i,"pvalint"]<=0.05){
marker <- as.character(QTLinfo[i,"topmarker"])
mdata<-SavedMarker[marker,]
WithVariant<- names(mdata)[mdata=="D"]
WithoutVariant<- names(mdata)[mdata=="B"]
Metabo<-gsub('_','.',QTLinfo[i,"Phenotype"])
main<-Metabo
plotFCMetabo(WithVariant,WithoutVariant,Metabo,main)
}
}